FOR · Laboratori

Grafi: cammino minimo, albero ricoprente e TSP

Lab 2: formulazioni di flusso, cut set e subtour elimination su grafi casuali
≈ 18 min di lettura · 3880 parole

I tre problemi di questo laboratorio sono i classici dell’ottimizzazione su grafi: il cammino minimo fra due nodi, l’albero ricoprente di costo minimo e il circuito hamiltoniano di costo minimo (il TSP). Il corso li affronta con algoritmi dedicati (Dijkstra, Prim, Kruskal); qui invece li risolviamo tutti con lo stesso strumento, una formulazione di programmazione lineare o lineare intera passata a Python-MIP. Il punto non è l’efficienza, che gli algoritmi specializzati vincono senza gara, ma imparare a scrivere vincoli indicizzati su nodi e archi e a gestire famiglie di vincoli esponenzialmente grandi.

I tre problemi sono legati: il cammino minimo ha una formulazione di flusso con variabili continue; l’albero ricoprente una formulazione con cut set e variabili binarie; il TSP una formulazione con vincoli di grado e subtour elimination, che si aggiungono iterativamente. Le istanze sono generate a caso, quindi ognuno ottiene un grafo diverso; in questa guida i risultati numerici si riferiscono a istanze con seme fissato, indicato ogni volta.

1. Generare un’istanza casuale#

Tutti e tre i notebook cominciano allo stesso modo: si generano kk punti con coordinate casuali nel quadrato [0,100]2[0, 100]^2, si calcolano le distanze euclidee e si collegano i punti abbastanza vicini.

import numpy as np
import math

k = 15
grid_size = 100          # lato del quadrato
d_max = .35 * grid_size  # oltre questa distanza non si crea l'arco

s = 0                    # sorgente
t = 12                   # destinazione

# np.random.seed(12345) # con un seme fissato l'istanza è riproducibile

point = grid_size * np.random.random((k, 2))
V = [i for i in range(k)]

np.random.random((k, 2)) produce una matrice k×2k \times 2 di numeri in [0,1)[0, 1); moltiplicandola per 100 si ottengono le coordinate. Il punto ii è point[i], un vettore di due componenti, e point[i, 0] e point[i, 1] sono le sue coordinate. Senza seed ogni esecuzione genera punti diversi: è normale che in aula ognuno vedesse un grafo diverso.

1.1 La matrice delle distanze#

La distanza fra tutti i punti si calcola in una riga:

distance = np.array([[math.sqrt(np.sum((point[i] - point[j])**2)) for i in V] for j in V])

Conviene smontarla dall’interno verso l’esterno, come fa il notebook stesso. Per due punti p1=(2,7)p_1 = (2, 7) e p2=(5,3)p_2 = (5, 3): p1 - p2 è il vettore differenza (−3,4)(-3, 4); (p1 - p2)**2 eleva al quadrato ogni componente, (9,16)(9, 16); np.sum somma, 2525; math.sqrt estrae la radice, 55. La doppia list comprehension [[... for i in V] for j in V] ripete il calcolo per ogni coppia, esattamente come [[i*j for i in range(1,6)] for j in range(1,6)] costruisce la tavola pitagorica; np.array trasforma la lista di liste in una matrice k×kk \times k, con distance[i, i] sempre zero.

Si può fare più in fretta?

Per migliaia di punti la doppia comprensione è lenta perché chiama math.sqrt k2k^2 volte da Python. Il notebook mostra la versione vettorizzata basata sul cosiddetto Euclidean trick: ∥pi−pj∥2=∥pi∥2+∥pj∥2−2 pi⋅pj\|p_i - p_j\|^2 = \|p_i\|^2 + \|p_j\|^2 - 2\, p_i \cdot p_j, che si calcola con un solo prodotto matriciale.

def slow_distance_matrix(points):
    return np.array([[math.sqrt(np.sum((pi - pj)**2)) for pi in points] for pj in points])

def fast_distance_matrix(x):
    xy = x @ x.T                       # tutti i prodotti scalari p_i · p_j
    x2 = xy.diagonal()[:, np.newaxis]  # le norme al quadrato, come colonna
    return np.abs(x2 + x2.T - 2. * xy)**0.5

# confronto (attenzione ai confronti fra numeri in virgola mobile)
eps = 1e-5
assert (np.abs(slow_distance_matrix(points) - fast_distance_matrix(points)) <= eps).all()

Con 1000 punti la versione lenta impiega secondi, quella veloce millisecondi. Per i lab la differenza è irrilevante.

1.2 Nodi, archi e disegno del grafo#

L’insieme degli archi contiene le coppie (i,j)(i, j) con i≠ji \ne j la cui distanza non supera dmax⁡d_{\max}:

import networkx as nx

A = [(i,j) for i in V for j in V if i != j and distance[i,j] <= d_max]

# coordinate come dizionario nodo -> (x, y), servono solo per il disegno
coord = {i: (point[i,0], point[i,1]) for i in V}

g = nx.Graph()
g.add_nodes_from(V)
g.add_edges_from(A)
nx.draw(g, pos=coord)

Due filtri, entrambi importanti. Senza i != j ogni nodo avrebbe un arco verso se stesso (distanza zero). Senza il filtro sulla distanza il grafo sarebbe completo e il disegno illeggibile; con dmax⁡d_{\max} al 35% del lato restano pochi archi per nodo, con il 75% se ne filtrano pochissimi. Il rovescio della medaglia: con dmax⁡d_{\max} troppo piccolo qualche nodo può restare isolato, e se è la sorgente o la destinazione il problema del cammino minimo diventa inammissibile. Vale la pena guardare il disegno prima di risolvere.

Un dettaglio di vocabolario che tornerà utile. Una coppia ordinata (i,j)(i, j) è un arco (arc), orientato da ii a jj; una coppia non ordinata {i,j}\{i, j\} è un lato (edge). Nell’insieme A qui sopra compaiono sia (i,j)(i, j) sia (j,i)(j, i), quindi il grafo è orientato con archi in entrambe le direzioni; nei problemi successivi si terrà solo la coppia con i<ji < j.

2. Il cammino minimo come problema di flusso#

Problema 1

Trovare il cammino minimo dal nodo 0 al nodo 12 su un grafo casuale con k=15k = 15 nodi, in cui due nodi sono collegati se la loro distanza euclidea non supera dmax⁡d_{\max}: generare i punti, creare gli archi, scrivere il modello di flusso, risolverlo e ricavare il cammino.

2.1 L’idea: un’unità di flusso da ss a tt#

Immaginiamo di spedire un’unità di flusso dalla sorgente ss alla destinazione tt lungo gli archi, pagando la distanza dell’arco per ogni unità che lo attraversa. Il flusso più economico percorre un cammino minimo. Per formalizzarlo si introduce una variabile fij≥0f_{ij} \ge 0 per ogni arco, la quantità di flusso sull’arco, e si impone la conservazione del flusso in ogni nodo. I nodi sono di tre tipi:

situvwzuscente − entrante = 1bs = 1uscente − entrante = 0bi = 0uscente − entrante = −1bt = −1flusso unitario da s a t: esce da s, attraversa i nodi intermedi, entra in t
Formulazione di flusso del cammino minimo

Dato il grafo G=(V,A)G = (V, A) con lunghezze dijd_{ij}, min⁡∑(i,j)∈Adijfijs.t.∑j∈V:(i,j)∈Afij−∑j∈V:(j,i)∈Afji=bi∀i∈Vfij≥0∀(i,j)∈A\begin{array}{lll} \min & \displaystyle\sum_{(i,j) \in A} d_{ij} f_{ij}\\ \text{s.t.} & \displaystyle\sum_{j \in V: (i,j) \in A} f_{ij} - \sum_{j \in V: (j,i) \in A} f_{ji} = b_i & \forall i\in V\\ & f_{ij} \ge 0 & \forall (i,j) \in A \end{array} dove bi=1b_i = 1 se i=si = s, bi=−1b_i = -1 se i=ti = t e bi=0b_i = 0 altrimenti.

Il primo termine del vincolo somma gli archi che partono da ii (flusso uscente), il secondo quelli che arrivano in ii (flusso entrante). Osserviamo che c’è un solo tipo di vincolo, uno per nodo, e che le variabili sono continue: non serve dichiararle binarie, perché per questa struttura di vincoli il solver restituisce comunque una soluzione a valori 0 e 1 (i vertici del poliedro dei flussi sono interi). Il valore dell’obiettivo è la lunghezza del cammino.

In parole semplici

Il vincolo di conservazione dice: «in ss parte un’unità, in tt ne arriva una, in mezzo nulla si perde e nulla si crea». Qualunque insieme di archi che soddisfa questa regola contiene un cammino da ss a tt; minimizzare la lunghezza totale sceglie il cammino più corto e scarta tutto il resto.

2.2 Il codice#

Le variabili sono indicizzate da coppie, quindi si usa un dizionario con chiavi tuple, non una lista:

import mip

m = mip.Model()

f = {(i,j): m.add_var() for (i,j) in A}   # limiti predefiniti 0 e +inf

Il parametro bib_i è un dizionario con tutti zeri, corretto poi nei due nodi speciali. Il vincolo di conservazione si scrive per ogni nodo con due xsum, la prima sugli archi in uscita e la seconda su quelli in ingresso; il confronto è == perché il vincolo è di uguaglianza:

b = {i: 0 for i in V}
b[s] = 1
b[t] = -1

for i in V:
    m.add_constr(mip.xsum(f[i,j] for j in V if (i,j) in A) -   # flusso uscente da i
                 mip.xsum(f[j,i] for j in V if (j,i) in A) ==  # flusso entrante in i
                 b[i])

Il filtro if (i,j) in A è essenziale: f[i,j] esiste solo per gli archi creati. Nella prima somma il nodo ii è fissato dal ciclo e jj scorre le destinazioni; nella seconda jj scorre le origini degli archi che entrano in ii, per questo la variabile è f[j,i]. Infine l’obiettivo, l’ottimizzazione e la lettura del cammino:

m.objective = mip.minimize(mip.xsum(distance[i,j] * f[i,j] for (i,j) in A))
m.optimize()

print(m.objective_value)              # lunghezza del cammino minimo
for (i,j) in A:
    if f[i,j].x > 0.5:
        print(i, "--->", j)

La soglia > 0.5 (o > 0.001, o > 0.9: qualunque valore fra 0 e 1) serve perché i valori restituiti dal solver sono numeri in virgola mobile: un arco nel cammino ha f.x uguale a 1 a meno di errori di arrotondamento, gli altri hanno 0.

Il notebook confeziona il disegno della soluzione in una funzione, che disegna prima l’intero grafo, poi solo gli archi con flusso positivo in rosso e più spessi, poi la radice in verde:

def draw_solution(V, A, x):
    g = nx.Graph()

    g.add_nodes_from(V)                              # tutto il grafo, senza evidenziare
    g.add_edges_from([(i,j) for (i,j) in A])
    nx.draw(g, pos=coord)

    g.clear()                                        # solo gli archi della soluzione
    g.add_edges_from([(i,j) for (i,j) in A if f[i,j].x > 0.001])
    nx.draw(g, pos=coord, width=4, edge_color='red')

    g.clear()                                        # la radice, in verde
    g.add_node(0)
    nx.draw(g, pos={0: coord[0]}, node_color='green')

draw_solution(V, A, f)

2.3 Un’istanza di riferimento#

Con np.random.seed(12345), k=15k = 15 e dmax⁡=35d_{\max} = 35 il grafo ha 62 archi (31 coppie di nodi vicini, ciascuna in entrambe le direzioni). Il cammino minimo da 0 a 12 è 0→3→10→120 \to 3 \to 10 \to 12, di lunghezza 87.02. Lo stesso valore si ottiene con l’algoritmo di Dijkstra, ed è la verifica migliore che la formulazione sia corretta.

Il cammino minimo sull’istanza di riferimento (seme 12345, k=15, d_{\\max}=35). In grigio gli archi del grafo, in arancione il cammino 0 \\to 3 \\to 10 \\to 12 di lunghezza 87.02; la sorgente è il nodo verde, la destinazione il nodo rosso.

3. L’albero ricoprente di costo minimo con i cut set#

Problema 2

Scrivere e risolvere la formulazione con cut set del problema dell’albero ricoprente di costo minimo (MST) su un grafo casuale.

3.1 Alberi ricoprenti e tagli#

Un albero ricoprente (spanning tree) di un grafo G=(V,E)G = (V, E) è un sottoinsieme di lati che collega tutti i nodi senza formare cicli. Un grafo ne ha in genere moltissimi; il problema chiede quello di costo minimo, dove il costo è la somma delle lunghezze dei lati scelti.

La formulazione si basa sulla proprietà che rende un insieme di lati «connesso»: comunque si separi l’insieme dei nodi in due parti, almeno un lato scelto deve attraversare la separazione.

Taglio (cut set)

Dato S⊂VS \subset V, il taglio δ(S)\delta(S) è l’insieme dei lati (i,j)∈E(i, j) \in E con esattamente un estremo in SS: i∈Si \in S e j∉Sj \notin S, oppure j∈Sj \in S e i∉Si \notin S.

S = {1, 2}1234567δ(S): (1,3) (1,4) (2,3) (2,7)altri lati

Se una soluzione non contiene nessun lato di δ(S)\delta(S), i nodi di SS sono scollegati dagli altri e non si tratta di un albero ricoprente. La condizione deve valere per ogni sottoinsieme SS non vuoto e diverso da VV: prendendo S={1}S = \{1\} si impone che il nodo 1 abbia almeno un lato, con S={1,2}S = \{1, 2\} che la coppia sia collegata al resto, e così via.

Formulazione con cut set dell'MST

Con una variabile binaria xijx_{ij} per ogni lato, uguale a 1 se il lato è nell’albero, min⁡∑(i,j)∈Ecijxijs.t.∑(i,j)∈δ(S)xij≥1∀S⊂V:S≠∅xij∈{0,1}∀(i,j)∈E\begin{array}{lll} \min & \displaystyle\sum_{(i,j)\in E} c_{ij} x_{ij}\\ \text{s.t.} & \displaystyle\sum_{(i,j)\in \delta(S)} x_{ij} \ge 1 & \forall S\subset V:S\neq \emptyset\\ & x_{ij} \in \{0, 1\} & \forall (i,j) \in E \end{array}

Rispetto al cammino minimo cambiano due cose: le variabili sono binarie (un lato c’è o non c’è, non è un flusso) e il grafo è non orientato, quindi si tiene una sola copia di ogni lato, con la convenzione i<ji < j. Osserviamo che nessun vincolo impone il numero di lati o l’assenza di cicli: con costi positivi, la minimizzazione elimina da sola ogni lato superfluo, e la soluzione ha esattamente k−1k - 1 lati.

Quanti vincoli sono? Uno per ogni sottoinsieme proprio non vuoto di VV. Ogni nodo può stare dentro o fuori da SS, quindi i sottoinsiemi sono 2k2^k; tolti l’insieme vuoto e VV stesso restano 2k−22^k - 2 vincoli. Con k=11k = 11 sono 2046; con k=30k = 30 sarebbero più di un miliardo. La formulazione è elegante, un solo tipo di vincolo, ma esponenziale: si può usare solo su grafi piccoli, ed è per questo che il notebook riduce kk a 11.

3.2 Il codice#

La generazione dell’istanza è la stessa, con k=11k = 11, un seme fissato e dmax⁡d_{\max} al 75% del lato, così che quasi tutti i lati siano presenti. L’insieme dei lati tiene solo i<ji < j:

import numpy as np
import math

k = 11
grid_size = 100
d_max = .75 * grid_size

np.random.seed(12345)

point = grid_size * np.random.random((k,2))
V = [i for i in range(k)]
distance = np.array([[math.sqrt(np.sum((point[i] - point[j])**2)) for i in V] for j in V])

# lati: i < j (una sola copia per coppia), e il filtro sulla distanza
E = [(i,j) for i in V for j in V if i < j and distance[i,j] <= d_max]

Tutti i sottoinsiemi di VV si generano con itertools: combinations(V, r) produce i sottoinsiemi di dimensione rr, chain.from_iterable li concatena per rr che varia in V, cioè da 0 a k−1k-1:

from itertools import chain, combinations
powerset = list(chain.from_iterable(combinations(V, r) for r in V))

Osserviamo che r in V fa variare la dimensione da 0 a k−1k - 1: l’insieme VV intero (dimensione kk) non viene mai generato, mentre l’insieme vuoto sì. Il modello scarta entrambi gli estremi con un if, e la condizione len(S) < k è in realtà già garantita; non fa male tenerla.

import mip

m = mip.Model()

x = {(i,j): m.add_var(var_type=mip.BINARY) for (i,j) in E}

for S in powerset:
    if len(S) > 0 and len(S) < k:
        m.add_constr(mip.xsum(x[i,j] for (i,j) in E
                              if (i in S and j not in S) or
                                 (j in S and i not in S)) >= 1)

m.objective = mip.minimize(mip.xsum(distance[i,j]*x[i,j] for (i,j) in E))

m.optimize()

sol = [(i,j) for (i,j) in E if x[i,j].x > 0.5]
print(f"{len(sol)} edges in the solution: {sol}")

La condizione dentro la xsum è la traduzione letterale di «(i,j)∈δ(S)(i, j) \in \delta(S)»: un estremo in SS e l’altro fuori, in uno dei due ordini. Il ciclo for S in powerset è il «∀S\forall S» della formulazione, ed è la riga che genera i 2046 vincoli.

3.3 Un’istanza di riferimento#

Con il seme 12345 e k=11k = 11 il grafo ha 43 lati. L’albero ottimo ha 10 lati, (0,3),(1,6),(1,7),(2,7),(2,10),(3,4),(3,9),(4,10),(5,8),(8,9)(0,3), (1,6), (1,7), (2,7), (2,10), (3,4), (3,9), (4,10), (5,8), (8,9), e costo 215.81; lo stesso albero si ottiene con l’algoritmo di Kruskal. Il disegno con draw_solution (adattato a x e alle coordinate point) mostra un insieme di lati che tocca ogni nodo e non chiude alcun ciclo.

L’albero ricoprente minimo dell’istanza di riferimento (seme 12345, k=11, d_{\\max}=75). In grigio i 43 lati del grafo, in arancione i 10 lati dell’albero ottimo, di costo totale 215.81.

4. Il TSP: vincoli di grado e subtour elimination#

Problema 3

Trovare il circuito hamiltoniano più corto attraverso nn nodi, note le distanze fra ogni coppia. Generare kk punti casuali, considerare il grafo completo, scrivere il modello con i soli vincoli di conservazione del flusso (senza subtour elimination), risolverlo e poi aggiungere iterativamente vincoli di subtour elimination finché la soluzione è un unico circuito.

4.1 La formulazione completa#

Un circuito hamiltoniano visita ogni nodo esattamente una volta e torna al punto di partenza. Con la variabile binaria xijx_{ij} per ogni lato (i<ji < j, grafo completo), ogni nodo deve essere toccato da esattamente due lati del circuito: uno per arrivare e uno per ripartire. È l’analogo, per un circuito, della conservazione del flusso della sezione 2, con «entrante più uscente uguale a 2» al posto di «uscente meno entrante uguale a bib_i».

Formulazione del TSP con subtour elimination

min⁡∑(i,j)∈Edijxijs.t.∑j∈V:(i,j)∈Exij=2∀i∈V∑(i,j)∈E:i∈C,j∈Cxij≤∣C∣−1∀C⊂V:C≠∅xij∈{0,1}∀(i,j)∈E\begin{array}{lll} \min & \displaystyle\sum_{(i,j) \in E} d_{ij} x_{ij}\\ \text{s.t.} & \displaystyle\sum_{j \in V: (i,j) \in E} x_{ij} = 2 & \forall i\in V\\ & \displaystyle\sum_{(i,j) \in E: i\in C, j\in C} x_{ij} \le |C| - 1 & \forall C\subset V: C \neq \emptyset\\ & x_{ij} \in \{0,1\} & \forall (i,j) \in E \end{array}

I vincoli di grado (la prima famiglia) da soli non bastano. Una soluzione che li rispetta può essere formata da più cicli disgiunti, i subtour: per esempio un triangolo su tre nodi e un quadrilatero sugli altri quattro, con ogni nodo di grado 2. Nulla nei vincoli di grado lo proibisce. La seconda famiglia lo esclude: per ogni sottoinsieme CC di nodi, i lati scelti interni a CC devono essere al più ∣C∣−1|C| - 1, perché un ciclo che passa per tutti e soli i nodi di CC ne userebbe esattamente ∣C∣|C|. Il ciclo su C={0,2,4}C = \{0, 2, 4\} usa i tre lati (0,2),(2,4),(0,4)(0,2), (2,4), (0,4); imporre x02+x24+x04≤2x_{02} + x_{24} + x_{04} \le 2 costringe almeno uno dei tre a uscire, e quindi il circuito a lasciare CC.

Anche qui i vincoli di subtour elimination sono esponenziali, uno per ogni sottoinsieme di nodi. Aggiungerli tutti in partenza è possibile solo per grafi minuscoli. La strategia del laboratorio è diversa e molto più generale.

4.2 Risolvere per iterazioni#

Si risolve il modello con i soli vincoli di grado. Se la soluzione è un unico circuito, è ottima anche per il problema completo (soddisfa vincoli che non abbiamo nemmeno scritto). Altrimenti si identificano i subtour nella soluzione, si aggiunge un vincolo di subtour elimination per ciascuno e si riottimizza; si ripete finché non restano subtour. In gergo i vincoli aggiunti a posteriori sono tagli e la procedura è una separazione di vincoli violati: è lo stesso schema che sta dietro ai risolutori industriali del TSP.

0123456solo vincoli di gradodue subtour, costo 210.150123456dopo due taglicircuito hamiltoniano, costo 224.04

Il codice della prima soluzione. L’istanza ha k=7k = 7 nodi e grafo completo, i lati sono le coppie con i<ji < j:

import numpy as np, math, networkx as nx, mip

k = 7
grid_size = 100
point = grid_size * np.random.random((k,2))
V = [i for i in range(k)]
distance = np.array([[math.sqrt(np.sum((point[i] - point[j])**2)) for i in V] for j in V])

E = [(i,j) for i in V for j in V if i < j]      # grafo completo, non orientato
coord = {i: (point[i,0], point[i,1]) for i in V}

m = mip.Model()
x = {(i,j): m.add_var(var_type=mip.BINARY) for (i,j) in E}

for i in V:
    m.add_constr(mip.xsum(x[j,i] for j in V if j < i) +
                 mip.xsum(x[i,j] for j in V if j > i) == 2)

m.objective = mip.minimize(mip.xsum(distance[i,j] * x[i,j] for (i,j) in E))
m.optimize()

Il vincolo di grado richiede attenzione agli indici: poiché esistono solo le variabili con primo indice minore del secondo, i lati incidenti in ii sono x[j,i] per j<ij < i e x[i,j] per j>ij > i. Scrivere x[i,j] per ogni jj darebbe un errore di chiave sul dizionario. Il disegno della soluzione riusa draw_solution con with_labels=True per mostrare i numeri dei nodi, e a questo punto in genere compaiono i subtour. Nell’istanza vista in aula erano il triangolo {0,2,4}\{0, 2, 4\} e il ciclo {6,3,1,5}\{6, 3, 1, 5\}; nel notebook l’esempio a mano è

m.add_constr(x[0,3] + x[3,4] + x[0,4] <= 2)
m.optimize()
draw_solution(V, E, x)

che spezza il triangolo {0,3,4}\{0, 3, 4\}. Nell’istanza di riferimento di questa guida (seme 6) la prima soluzione ha i subtour {0,1,3}\{0, 1, 3\} e {2,4,5,6}\{2, 4, 5, 6\} con costo 210.15; aggiunto il taglio su {0,1,3}\{0, 1, 3\} la nuova soluzione costa 212.65 e ha i subtour {0,1,3,6}\{0, 1, 3, 6\} e {2,4,5}\{2, 4, 5\}; aggiunto il taglio su {0,1,3,6}\{0, 1, 3, 6\} si ottiene finalmente il circuito 0,1,3,2,4,5,6,00, 1, 3, 2, 4, 5, 6, 0 di costo 224.04. Tre risoluzioni al posto dei 27−2=1262^7 - 2 = 126 vincoli della formulazione completa.

4.3 Automatizzare: trovare i cicli#

Aggiungere i tagli a mano è istruttivo ma non praticabile. Serve una funzione che, data la soluzione, restituisca i cicli. Il notebook costruisce la lista di adiacenza dei lati scelti e la visita in profondità (DFS) colorando i nodi: 0 non visitato, 1 in visita, 2 visitato; quando la visita incontra un nodo già in stato 1 ha trovato un ciclo e lo ricostruisce risalendo i padri.

def get_cycles(V, E, x):
    graph = [[] for i in V]

    for (i,j) in E:
        if x[i,j].x > 0.5:
            graph[i].append(j)
            graph[j].append(i)

    cycles = []
    color = [0 for i in V]
    par = [-1 for i in V]

    for i in V:
        if par[i] == -1:
            dfs_cycle(graph, cycles, i, -1, color, par)

    return [c for c in cycles if len(c) < len(V) and len(c) > 0]

def dfs_cycle(graph, cycles, u, p, color, par):
    if color[u] == 2:            # nodo già completamente visitato
        return

    if color[u] == 1:            # nodo in visita: ciclo trovato, lo ricostruisco dai padri
        v = []
        cur = p
        v.append(cur)
        while cur != u:
            cur = par[cur]
            v.append(cur)
        cycles.append(v)
        return

    par[u] = p
    color[u] = 1                 # visita in corso

    for v in graph[u]:
        if v == par[u]:          # non tornare indietro sul lato appena percorso
            continue
        dfs_cycle(graph, cycles, v, u, color, par)

    color[u] = 2                 # visita completata

get_cycles restituisce solo i cicli con meno di ∣V∣|V| nodi: un ciclo che li contiene tutti è un circuito hamiltoniano, e in quel caso la lista è vuota. Il ciclo principale risolve, cerca i subtour, aggiunge un taglio per il primo trovato e ripete:

m = mip.Model()

x = {(i,j): m.add_var(var_type=mip.BINARY, name="x{}{}".format(str(i), str(j))) for (i,j) in E}

for i in V:
    m.add_constr(mip.xsum(x[j,i] for j in V if j < i) +
                 mip.xsum(x[i,j] for j in V if j > i) == 2)

m.objective = mip.minimize(mip.xsum(distance[i,j] * x[i,j] for (i,j) in E))
m.optimize()

cycles = get_cycles(V, E, x)
while len(cycles) > 0:
    cycle = cycles.pop()
    cycle_edges = [x[i,j] for (i,j) in E if i in cycle and j in cycle]
    c = m.add_constr(mip.xsum(cycle_edges) <= len(cycle) - 1)
    print("Added constr:")
    print(c)
    m.optimize()
    cycles = get_cycles(V, E, x)

print(m.objective_value)
draw_solution(V, E, x)

Il name= dato alle variabili serve solo a rendere leggibile la stampa del vincolo aggiunto, per esempio x03 + x34 + x04 <= 2.

Una correzione al notebook della soluzione

Nel notebook il termine di destra del taglio è len(cycle_edges) - 1, dove cycle_edges sono tutti i lati del grafo completo con entrambi gli estremi nel ciclo, cioè (∣C∣2)\binom{|C|}{2} lati. Per un ciclo di 3 nodi i lati interni sono 3 e il vincolo è ≤2\le 2, corretto. Per un ciclo di 4 nodi i lati interni sono 6 e il vincolo diventa ≤5\le 5: il subtour ne usa 4 e lo soddisfa, quindi il taglio non lo elimina e il ciclo while può ripetersi all’infinito sulla stessa soluzione. Il termine di destra corretto, quello della formulazione, è ∣C∣−1|C| - 1, ovvero len(cycle) - 1, come scritto nel listato qui sopra.

La ricerca dei cicli non ha bisogno del solver: dato l’insieme dei lati scelti, si può provare direttamente nella pagina. Il blocco seguente usa una variante più semplice, che raccoglie le componenti connesse del grafo dei lati scelti (con i vincoli di grado ogni componente è un ciclo).

In parole semplici

Il modello «rilassato» con i soli vincoli di grado è ottimista: trova soluzioni a pezzi che costano meno del vero ottimo. Ogni taglio gli vieta uno specifico modo di andare a pezzi, il costo sale un poco, e dopo qualche giro l’unica soluzione rimasta è un circuito unico. Non si scrivono mai tutti i vincoli, solo quelli che la soluzione corrente viola.

5. Riepilogo e glossario#

Problema Variabili Vincoli Numero di vincoli Trucco
Cammino minimo fij≥0f_{ij} \ge 0 continue, su archi orientati conservazione del flusso, uno per nodo ∣V∣\lvert V \rvert LP basta: la soluzione è già intera
MST xij∈{0,1}x_{ij} \in \{0,1\} su lati, i<ji<j almeno un lato per ogni taglio δ(S)\delta(S) 2k−22^k - 2 si generano con itertools, solo per kk piccolo
TSP xij∈{0,1}x_{ij} \in \{0,1\} su lati, grafo completo grado 2 in ogni nodo; subtour elimination ∣V∣\lvert V \rvert + esponenziali i tagli si aggiungono solo quando violati

In tutti e tre i casi il codice ha la stessa forma: dati generati con NumPy, variabili in un dizionario indicizzato da coppie, un ciclo for per la famiglia di vincoli, xsum per le somme, x[i,j].x > 0.5 per leggere la soluzione.

Termine Significato
arco / lato coppia ordinata (i,j)(i,j) (orientato) / coppia non ordinata {i,j}\{i,j\}
conservazione del flusso in ogni nodo, flusso uscente meno entrante uguale a bib_i
cut set δ(S)\delta(S) lati con esattamente un estremo in SS
albero ricoprente sottoinsieme di lati connesso e senza cicli che tocca tutti i nodi
circuito hamiltoniano ciclo che visita ogni nodo esattamente una volta
subtour ciclo che visita solo una parte dei nodi
taglio (cut) vincolo aggiunto al modello perché violato dalla soluzione corrente
separazione procedura che trova un vincolo violato (qui, la ricerca dei cicli)

Elaborato dai tre notebook del Lab 2 (shortest path, minimum spanning tree cutset, tsp) e dalla registrazione della seconda lezione di laboratorio.